/**
 * 
 */
package simulation;

//import static common.FSave.fSave;

import org.apache.commons.math.stat.StatUtils;

import java.util.ArrayList;
import java.util.Date;

/**
 * @author vladimir
 * 
 */
public class Main {

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// Parameters
		// Neighbors in spatial direction
		int[][] neighbors = { { 0, 1 }, { 0, -1 }, { 1, 0 }, { -1, 0 },
				{ -1, -1 }, { 1, 1 } };

		String modelName = "triang";
		// int[][] neighbors = { { 0, 1 }, { 0, -1 }, { 1, 0 }, { -1, 0 } };
		// String modelName = "square";

		// int[][] neighbors = { { 1 }, { -1 } };
		// String modelName = "chain";

		// spatial coupling constant
		double J = -1;

		// Array
		// describing
		// size
		// of
		// the
		// lattice
		int[] lRange = { 12, 15, 18, 24, 30, 36, 48, 60, 75, 90};

		// Array describing temperature
		double[] tRange = common.Range.range(0.05, 1.0, 0.05);
		// double[] tRange = { 0.05 };
		// Array describing Magnetic field
		// double[] bRange = common.Range.range(0.0, 5.0, 0.1);
//		double[] bRange = { 0.1, 0.2, 0.4, 0.8, 1.2, 1.6};
        double[] bRange = { 0.4};

		int dimension = neighbors[0].length;

		boolean info = true;

		for (int L : lRange) {
			for (double B : bRange) {// go through different magnetic field
										// values
				double tB = Math.round(1000 * B) / 1000.0d;
				for (double T : tRange) {// go through different temperatures
					double tT = Math.round(1000 * T) / 1000.0d;
					java.util.Date Idate = new Date();

					CorrTime Cor = new CorrTime(neighbors, J, L, B, T,
							modelName);

					Lattice Lat = Cor.Lat;

					int MCS = -1;
					int Thermal_MCS = -1;
					if (Cor.ThermalMCS() > Cor.mcs) {
						// number of the Monte Carlo Sweeps
						MCS = Cor.MCS() - Cor.mcs;
						// number of the Monte Carlo Sweeps to reach equilibrium
						Thermal_MCS = Cor.ThermalMCS() - Cor.mcs;
					} else {
						// number of the Monte Carlo Sweeps
						MCS = Cor.MCS() - Cor.ThermalMCS();
						// number of the Monte Carlo Sweeps to reach equilibrium
						Thermal_MCS = 0;
					}

					// Correlation time
					int sweepStep = Cor.SweepStep();

					ArrayList<Double> tE = new ArrayList<Double>();
					ArrayList<Double> tM = new ArrayList<Double>();
					ArrayList<Double> tM2 = new ArrayList<Double>();
					ArrayList<Double> tM4 = new ArrayList<Double>();
                    ArrayList<Double> tM_total = new ArrayList<Double>();
					ArrayList<Double> tM2_total = new ArrayList<Double>();
					ArrayList<Double> tM4_total = new ArrayList<Double>();

					// This arrays to keep stuff for triangular
					// antiferromagnetic
					ArrayList<Double> tOP = new ArrayList<Double>();
					ArrayList<Double> tOPT = new ArrayList<Double>();

					ArrayList<Double> tOP2 = new ArrayList<Double>();
					ArrayList<Double> tOPT2 = new ArrayList<Double>();

					ArrayList<Double> tOP4 = new ArrayList<Double>();
					ArrayList<Double> tOPT4 = new ArrayList<Double>();

					ArrayList<Double> tOP6 = new ArrayList<Double>();
					ArrayList<Double> tOPT6 = new ArrayList<Double>();

					ArrayList<Double> tOP6cos = new ArrayList<Double>();
					ArrayList<Double> tOPT6cos = new ArrayList<Double>();

					ArrayList<Double> nCuts = new ArrayList<Double>();

					if (info) {
						System.out.println();
						System.out.println("modelName : " + modelName);
						System.out.println("J - " + J);
						System.out.println("L - " + L);
						System.out.println("B - " + tB);
						System.out.println("T - " + tT);
						System.out.println("Monte Carlo Sweeps - " + Cor.MCS());

						System.out.println("Thermal Monte Carlo Sweeps - "
								+ Cor.ThermalMCS());
						System.out.println("Correlation time - " + sweepStep);
					}
					//path to save data
					String resultData = common.PathJoin.pathJoin("..","results", modelName, "J_"
						+ J, "B_" + tB, "L_" + L);

					common.FSave.fSave(common.PathJoin.pathJoin(resultData,	"autoCorrelationTime"), 
						"T" + tT + "_" + Idate.getTime(), Cor.SweepStep());

					for (int sweep = 0; sweep < MCS; sweep++) {
						Lat.update();
						if ((Thermal_MCS < sweep) && (sweep % sweepStep == 0)) {
							tE.add(Lat.Energy());

							nCuts.add(Lat.nCuts());

							tM.add(Lat.M(1));
							tM2.add(Lat.M(2));
							tM4.add(Lat.M(4));
							double M_total = Lat.M_total();
							tM_total.add(Math.pow(M_total, 1));
							tM2_total.add(Math.pow(M_total, 2));
							tM4_total.add(Math.pow(M_total, 4));

							double oP = Lat.orderParameter();
							double oPT = Lat.orderParameter_total();

							tOP.add(oP);
							tOPT.add(oPT);

							tOP2.add(Math.pow(oP, 2));
							tOPT2.add(Math.pow(oPT, 2));

							tOP4.add(Math.pow(oP, 4));
							tOPT4.add(Math.pow(oPT, 4));

							tOP6.add(Math.pow(oP, 6));
							tOPT6.add(Math.pow(oPT, 6));

							tOP6cos.add(Math.pow(oP, 6) * Lat.cos6teta());
							tOPT6cos.add(Math.pow(oPT, 6) * Lat.cos6teta_total());

						}
					}

					// save data to file
                    
					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "Energy"), "T"
									+ tT + "_" + Idate.getTime(),
							StatUtils.mean(common.Todouble.todouble(tE)));

					double M2 = StatUtils.mean(common.Todouble.todouble(tM2));

					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "M2"), "T" + tT
									+ "_" + Idate.getTime(), M2 );

					double M4 = StatUtils.mean(common.Todouble.todouble(tM4));

					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "M4"), "T" + tT
									+ "_" + Idate.getTime(), M4 );

					double U_L = 1 - M4 / (3 * Math.pow(M2, 2));

					common.FSave.fSave(common.PathJoin.pathJoin(resultData, "U_L"), "T" + tT + "" + Idate.getTime(), U_L);

					double X = Math.pow(L, dimension) * StatUtils.variance(common.Todouble.todouble(tM)) / T;

					common.FSave.fSave(common.PathJoin.pathJoin(resultData, "X"), "T" + tT + "_" + Idate.getTime(), X );


					double M2_total = Math
							.abs(StatUtils.mean(common.Todouble
									.todouble(tM2_total)));

					common.FSave.fSave(common.PathJoin.pathJoin(resultData,
							"M2_total"), "T" + tT + "" + Idate.getTime(), M2_total);

					double M4_total = StatUtils
											.mean(common.Todouble
													.todouble(tM4_total));
					common.FSave
							.fSave(common.PathJoin.pathJoin(resultData,
									"M4_total"),
									"T" + tT + "_" + Idate.getTime(), M4_total);

					double U_Ltotal = 1 - M4_total / (3 * Math.pow(M2_total, 2));						

					common.FSave.fSave(common.PathJoin.pathJoin(resultData, "U_Ltotal"), "T" + tT + "" + Idate.getTime(), U_Ltotal);

					double X_total = Math.pow(L, dimension) * StatUtils.variance(common.Todouble.todouble(tM_total)) / T;

					common.FSave.fSave(common.PathJoin.pathJoin(resultData, "X_total"), "T" + tT + "_" + Idate.getTime(), X_total );

//					common.FSave.fSave(
//							common.PathJoin.pathJoin(resultData, "OP"), "T" + tT
//									+ "_" + Idate.getTime(),
//							StatUtils.mean(common.Todouble.todouble(OP)));

                    // common.FSave.fSave(
                    //         common.PathJoin.pathJoin(resultData, "X_OP"),
                    //         "T" + tT + "_" + Idate.getTime(),
                    //         Math.pow(L, dimension)
                    //                 * StatUtils.variance(common.Todouble
                    //                 .todouble(OP)) / T);


					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OPT"), "T" + tT
									+ "_" + Idate.getTime(),
							StatUtils.mean(common.Todouble.todouble(tOPT)));

					double OP2 = StatUtils.mean(common.Todouble.todouble(tOP2));
					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OP2"), "T" + tT
									+ "_" + Idate.getTime(), OP2 );
					double OP4 = StatUtils.mean(common.Todouble.todouble(tOP4));
					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OP4"), "T" + tT
									+ "_" + Idate.getTime(), OP4 );

					double U_LOP = 1 - OP4 / (3 * Math.pow(OP2, 2));

					common.FSave.fSave(common.PathJoin.pathJoin(resultData, "U_LOP"), "T" + tT + "" + Idate.getTime(), U_LOP);

					double X_OP = Math.pow(L, dimension)
                                    * StatUtils.variance(common.Todouble
                                    .todouble(tOP)) / T;
					common.FSave.fSave(
                            common.PathJoin.pathJoin(resultData, "X_OP"),
                            "T" + tT + "_" + Idate.getTime(), X_OP );


                    double V_LOP = Math.pow(L, dimension) * OP2 / X_OP;

                    common.FSave.fSave(common.PathJoin.pathJoin(resultData,
                            "V_LOP"),
                            "T" + tT + "_" + Idate.getTime(), V_LOP);

					double OPT2 = StatUtils.mean(common.Todouble.todouble(tOPT2));
					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OPT2"), "T" + tT
									+ "_" + Idate.getTime(), OPT2 );


					double OPT4 = StatUtils.mean(common.Todouble.todouble(tOPT4));
					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OPT4"), "T" + tT
									+ "_" + Idate.getTime(), OPT4);

					double U_LOPT = 1 - OPT4 / (3 * Math.pow(OPT2, 2));

					common.FSave.fSave(common.PathJoin.pathJoin(resultData, "U_LOPT"), "T" + tT + "" + Idate.getTime(), U_LOPT);

					double X_OPT = Math.pow(L, dimension)
                                    * StatUtils.variance(common.Todouble
                                    .todouble(tOPT)) / T;
					common.FSave.fSave(
                            common.PathJoin.pathJoin(resultData, "X_OPT"),
                            "T" + tT + "_" + Idate.getTime(), X_OPT );

					double V_LOPT = Math.pow(L, dimension) * OPT2 / X_OPT;

                    common.FSave.fSave(common.PathJoin.pathJoin(resultData,
                            "V_LOPT"),
                            "T" + tT + "_" + Idate.getTime(), V_LOPT);

					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OP6"), "T" + tT
									+ "_" + Idate.getTime(),
							StatUtils.mean(common.Todouble.todouble(tOP6)));
					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OPT6"), "T" + tT
									+ "_" + Idate.getTime(),
							StatUtils.mean(common.Todouble.todouble(tOPT6)));
					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OP6cos"), "T"
									+ tT + "_" + Idate.getTime(),
							StatUtils.mean(common.Todouble.todouble(tOP6cos)));
					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "OPT6cos"), "T"
									+ tT + "_" + Idate.getTime(),
							StatUtils.mean(common.Todouble.todouble(tOPT6cos)));
					common.FSave.fSave(common.PathJoin.pathJoin(resultData,
							"autoCorrelationTime"),
							"T" + tT + "_" + Idate.getTime(), Cor.SweepStep());

					double C = StatUtils.variance(common.Todouble.todouble(tE))
							/ Math.pow(T, 2) * Math.pow(L, dimension)
							- StatUtils.mean(common.Todouble.todouble(nCuts));

					common.FSave.fSave(
							common.PathJoin.pathJoin(resultData, "C"), "T" + tT
									+ "_" + Idate.getTime(), C);

					if (info) {
						java.util.Date Fdate = new Date();
						long timeDiff = Fdate.getTime() - Idate.getTime();
						System.out.println("Waste of time: " + timeDiff);
					}
				}
			}
		}
	}
}
